Integration of multimodal cues does not alter mean but reduces variance in avian responses to predators: a systematic review and meta-analysis

Author

Kim + Shinichi et al.

Published

August 20, 2023

1 Setting up

Code
# install.packages("pacman")
# pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)

# devtools::install_github("daniel1noble/orchaRd", force = TRUE)
library(pacman)

p_load(
    tidyverse,
    here,
    lme4,
    orchaRd,
    metafor,
    patchwork,
    alluvial,
    ggalluvial,
    ape,
    easyalluvial,
    clubSandwich,
    emmeans,
    MuMIn,
    kableExtra,
    pander
)

# library(tidyverse)
# library(here)
# library(lme4)
# library(orchaRd)
# library(metafor)
# library(patchwork)
# library(alluvial)
# library(ggalluvial)
# library(easyalluvial)
# library(ape)
# library(clubSandwich)
# library(emmeans)
# library(MuMIn)
# library(kableExtra)
# library(pander)
# making metafor talk to MuMIn
eval(metafor:::.MuMIn)
# install.packages("pak")
# pak::pak("MichelNivard/gptstudio")

2 Getting data loaded

Code
#dat_full <- read.csv(here("data/dat_04_04_2023.csv"))
#dat_full <- read.csv(here("data/dat_28_06_2023.csv"))
dat_full <- read.csv(here("data/dat_19_07_2023.csv"))

# add phylogenetic tree - only topologies
# TODO? - we could get better tree from birdtree.org
# we can do 50 different trees as in 
# https://academic.oup.com/sysbio/article/68/4/632/5267840
tree_top <- read.tree(here("R/birds_MA.tre"))

# tree with branch lengths
tree <- compute.brlen(tree_top)
#plot(tree)
# turning it into a correlation matrix
cor_tree <- vcv(tree,corr=T)

3 Custom functions

Code
# custom functions

#' Title: Contrast name generator
#'
#' @param name: a vector of character strings
cont_gen <- function(name) {
  combination <- combn(name, 2)
  name_dat <- t(combination)
  names <- paste(name_dat[, 1], name_dat[, 2], sep = "-")
  return(names)
}

#' @title get_pred1: intercept-less model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred1 <- function (model, mod = " ") {
  name <- firstup(as.character(stringr::str_replace(row.names(model$beta), mod, "")))
  len <- length(name)
  
   if (len != 1) {
        newdata <- diag(len)
        pred <- metafor::predict.rma(model, 
                                     newmods = newdata,
                                     tau2.levels = 1:len)
    }
    else {
        pred <- metafor::predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title get_pred2: normal model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred2 <- function (model, mod = " ") {
  name <- as.factor(str_replace(row.names(model$beta), 
                                paste0("relevel", "\\(", mod,", ref = name",
                                       "\\)"),""))
  len <- length(name)
  
  if(len != 1){
  newdata <- diag(len)
  pred <- predict.rma(model, intercept = FALSE, newmods = newdata[,-1])
  }
  else {
    pred <- predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title mr_results
#' @description Function to put results of meta-regression and its contrasts
#' @param res1: data frame 1
#' @param res1: data frame 2
mr_results <- function(res1, res2) {
  restuls <-tibble(
    `Fixed effect` = c(as.character(res1$name), cont_gen(res1$name)),
    Estimate = c(res1$estimate, res2$estimate),
    `Lower CI [0.025]` = c(res1$lowerCL, res2$lowerCL),
    `Upper CI  [0.975]` = c(res1$upperCL, res2$upperCL),
    `P value` = c(res1$pval, res2$pval),
    `Lower PI [0.025]` = c(res1$lowerPR, res2$lowerPR),
    `Upper PI  [0.975]` = c(res1$upperPR, res2$upperPR),
  )
}


#' @title all_models
#' @description Function to take all possible models and get their results
#' @param model: intercept-less model
#' @param mod: the name of a moderator 

all_models <- function(model, mod = " ", type = "homo") {
  
  # getting the level names out
  level_names <- levels(factor(model$data[[mod]]))
  dat2 <- model$data
  mod <- mod

  VCV1 <- vcalc(vi = dat2$Vd,
             cluster = dat2$SubjectID,
             rho = 0.5)
  #model$data[[mod]] <- factor(model$data[[mod]], ordered = FALSE)
  # meta-regression: contrasts 
  # helper function to run metafor meta-regression
  run_rma1 <- function(name) {
    rma.mv(yi = SMD, 
         V = VCV1, 
         mods = ~ relevel(dat2[[mod]], ref = name), 
         random = list(~1 | FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat2)
   }

    run_rma2 <- function(name) {
    rma.mv(yi = SMD,
         V = VCV1,
         mods = ~ relevel(dat2[[mod]], ref = name),
         random = list(~1 | FocalSpL,
                             ~1 | RecNo,
                             formula(paste("~", mod, "| Obs_ID"))),
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat2)
   }

  # results of meta-regression including all contrast results; taking the last level out ([-length(level_names)])
  # this does not work for hetero model?
  if (type == "homo"){
  model_all <- map(level_names[-length(level_names)], run_rma1)
  } else {
  model_all <- map(level_names[-length(level_names)], run_rma2)
  }
  
  # getting estimates from intercept-less models (means for all the groups)
  res1 <- get_pred1(model, mod = mod)
  
  # getting estiamtes from all contrast models
  res2_pre <- map(model_all, ~ get_pred2(.x, mod = mod))
  
  # a list of the numbers to take out unnecessary contrasts
  contra_list <- Map(seq, from=1, to=1:(length(level_names) - 1))
  res2 <- map2_dfr(res2_pre, contra_list, ~.x[-(.y), ]) 
  # creating a table
  res_tab <- mr_results(res1, res2) %>% 
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")
  
  # results
  res_tab

}

# test 
#all_models(mod1, mod = "Treat_mod")
#all_models(mod1s, mod = "Treat_mod")
#all_models(mod1s2, mod = "Treat_mod")

4 Preparing data set

Code
# function for calculating variance
Vd_func <- function(d, n1, n2, design, r = 0.5){
  # independent design
  if(design == "among"){
    var <- (n1 + n2) / (n1*n2) + d^2 / (2 * (n1 + n2 - 2)) # variance
  } else { # dependent design
    var <- 2*(1-r) / n1 + d^2 / (2*(n1 - 1)) # variance
  }
  var # return variance
}

# getting Hedges' g - get small size bias corrected effect size
dat_full$SMD <- dat_full$d / (1 - 3/(4 * (dat_full$NTreat + dat_full$Ncontrol) - 9))

# flipping d 
dat_full$SMD <- dat_full$d*dat_full$Direction*dat_full$PredictedDirection


# calucating Vd
dat_full$Vd <- with(dat_full, pmap_dbl(list(SMD, NTreat, Ncontrol, Design), Vd_func))


# extra useful function
# function for getting mean and sd from median, quartiles and sample size
# get_mean_sd <- function(median, q1, q3, n){
#   sd <- (q3 - q1) / (2 * (qnorm((0.75 * n - 0.125) / (n + 0.25)))) # sd
#   mean <- (median + q1 + q3)/3 # mean
#   c(mean, sd)
# }


# observation id
dat_full$Obs_ID <- 1:nrow(dat_full)
dat_full$Phylo <- gsub(" ", "_", dat_full$FocalSpL)

# filtering very large variance and also very small sample size
dat_int <- dat_full %>% filter(Vd < 10 & Ncontrol > 2 & NTreat > 2)

#dim(dat_full)
#dim(dat_int)


# sorting out modality stuff
# creat - 1,2,3 modality - also easier classification A, O, V (AOV = L) 

dat_int %>% mutate(Treat_mod = case_when(Treatment == "A" ~ "A",
                                          Treatment == "AV" ~ "AV",
                                          Treatment == "AVG" ~ "AV",
                                          Treatment == "AVM" ~ "AV",
                                          Treatment == "L" ~ "AVO",
                                          Treatment == "O" ~ "O",
                                          Treatment == "OV" ~ "OV",
                                          Treatment == "V" ~ "V",
                                          Treatment == "VG" ~ "V",
                                          Treatment == "VM" ~ "V",
                                          Treatment == "VP" ~ "V"),
                    # into how many
                    Treat_No = case_when(Treatment == "A" ~ 1,
                                         Treatment == "AV" ~ 2,
                                         Treatment == "AVG" ~ 2,
                                         Treatment == "AVM" ~ 2,
                                         Treatment == "L" ~ 3,
                                         Treatment == "O" ~ 1,
                                         Treatment == "OV" ~ 2,
                                         Treatment == "V" ~ 1,
                                         Treatment == "VG" ~ 1,
                                         Treatment == "VM" ~ 1,
                                         Treatment == "VP" ~ 1),
                    # des it have some add-ons
                    Add_on = case_when(Treatment == "A" ~ "No",
                                         Treatment == "AV" ~ "No",
                                         Treatment == "AVG" ~ "Yes",
                                         Treatment == "AVM" ~ "Yes",
                                         Treatment == "L" ~ "No",
                                         Treatment == "O" ~ "No",
                                         Treatment == "OV" ~ "No",
                                         Treatment == "V" ~ "No",
                                         Treatment == "VG" ~ "Yes",
                                         Treatment == "VM" ~ "Yes",
                                         Treatment == "VP" ~ "Yes"),

                      ) -> dat

# creating data just for A, V, and AV 
dat_short <- dat %>% filter(Treat_mod == "A" | Treat_mod == "V" | Treat_mod == "AV")

# for add-on, we only need V and AV
dat_short_add <- dat %>% filter(Treat_mod == "AV" | Treat_mod == "V")


dat <- dat %>%
  mutate_if(is.character, as.factor)

# reordering factors for better visualization

# dat$Treat_mod <- factor(dat$Treat_mod, levels = c("A", "V", "AV", "O", "OV", "AVO"))

5 Exploratory visualization

For Treat_mod (Treatment), we will only visualise A, V, and AV as O $r $, OV, and AVO are much rarer. But for Type (Trait type), we will use all data.

Code
# reordering dat_shrot$Treat_mod for better visualisation
dat_short$Treat_mod <- factor(dat_short$Treat_mod, levels = c("A", "V", "AV"))

# Treatment vs Type
dat_short %>% group_by(Treat_mod, Type) %>%
  summarise(n = n()) -> tab

#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = Type)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and trait type")

Code
# tuning Treatment duration into a binary variable
dat_short %>% mutate(TDration = case_when(duration_days < 1 ~ "< 1 day",
                                          duration_days >= 1 ~ "> 1 day")) -> dat_short

# reformatting data
dat_short %>% group_by(Treat_mod, TDration) %>%
  summarise(n = n()) -> tab

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = TDration)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and duration of treatment")

Code
 # Treat_mod vs Design

dat_short %>% group_by(Treat_mod, Sex) %>%
  summarise(n = n()) -> tab
#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = Sex)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and study setting")

Code
# Treat_mod vs Design
dat_short %>% filter(PredTo != "") %>%
  group_by(Treat_mod, PredTo) %>%
  summarise(n = n()) -> tab

tab$PredTo <- factor(tab$PredTo, levels = c("A", "N", "B"), labels = c("Adult", "Nestling", "Both"))

#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = PredTo)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and study setting")

Code
 # Treat_mod vs Design

dat_short %>% group_by(Treat_mod, Design) %>%
  summarise(n = n()) -> tab
#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = Design)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and study setting")

Code
 # Treat_mod vs Design

dat_short %>% group_by(Treat_mod, Season) %>%
  summarise(n = n()) -> tab
#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = Season)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and study setting")

Code
 # Treat_mod vs Design

dat_short %>% group_by(Treat_mod, Setting) %>%
  summarise(n = n()) -> tab
#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = Setting)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and study setting")

Code
# Treat_mod vs ControlType

dat_short %>% filter(ControlType != "mix") %>% group_by(Treat_mod, ControlType) %>%
  summarise(n = n()) -> tab
#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = ControlType)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 6, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and control type")

Code
# Type vs duration_days

# turn duration_days into a binary factor: less than 1 and more than 1

dat %>% mutate(TDration = case_when(duration_days < 1 ~ "< 1 day",
                                          duration_days >= 1 ~ "> 1 day")) -> dat


dat %>% group_by(Type, TDration) %>%
  summarise(n = n()) -> tab2

# using ggaruvial
ggplot(tab2,
       aes(y = n,
           axis1 = Type,
           axis2 = TDration)) +
  geom_alluvium(aes(fill = Type)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Trait type and duration of treatment")

Code
# Type vs duration_days

# turn duration_days into a binary factor: less than 1 and more than 1

dat %>% mutate(TDration = case_when(duration_days < 1 ~ "< 1 day",
                                          duration_days >= 1 ~ "> 1 day")) -> dat

dat %>% group_by(Type, Sex) %>%
  summarise(n = n()) -> tab

# reordering
tab$Sex <- factor(tab$Sex, levels = c("F", "M", "both"), labels = c("Female", "Male", "Both"))


# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Type,
           axis2 = Sex)) +
  geom_alluvium(aes(fill = Type)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Trait type and duration of treatment")

6 Meta-analysis

Code
##| warning: false
# VCV matrix
# TODO - new phylogeny

VCV <- vcalc(vi = dat$Vd,
             cluster = dat$SubjectID,
             rho = 0.5)

mod_ma <- rma.mv(yi = SMD,
       V = VCV, 
       random = list(#~1 | Phylo,
                     ~1 | FocalSpL,
                     ~1 | RecNo,
                     ~1 | SubjectID, # incoprated as VCV
                     ~1 | Obs_ID),
      # R = list(Phylo = cor_tree),
       test = "t",
       method = "REML", 
       sparse = TRUE,
       data = dat)

summary(mod_ma)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-894.6453  1789.2906  1799.2906  1821.5901  1799.3853   

Variance Components:

            estim    sqrt  nlvls  fixed     factor 
sigma^2.1  0.0111  0.1053     90     no   FocalSpL 
sigma^2.2  0.1450  0.3808    116     no      RecNo 
sigma^2.3  0.0000  0.0001    153     no  SubjectID 
sigma^2.4  0.6872  0.8290    640     no     Obs_ID 

Model Results:

estimate      se    tval   df    pval   ci.lb   ci.ub      
  0.4139  0.0653  6.3401  639  <.0001  0.2857  0.5421  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# TODO - think about whether we add this or not
robust(mod_ma, cluster = dat$SubjectID)

Multivariate Meta-Analysis Model (k = 640; method: REML)

Variance Components:

            estim    sqrt  nlvls  fixed     factor 
sigma^2.1  0.0111  0.1053     90     no   FocalSpL 
sigma^2.2  0.1450  0.3808    116     no      RecNo 
sigma^2.3  0.0000  0.0001    153     no  SubjectID 
sigma^2.4  0.6872  0.8290    640     no     Obs_ID 

Number of estimates:   640
Number of clusters:    153
Estimates per cluster: 1-29 (mean: 4.18, median: 2)

Model Results:

estimate      se¹    tval¹   df¹    pval¹   ci.lb¹   ci.ub¹      
  0.4139  0.0540   7.6691   152   <.0001   0.3073   0.5206   *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1) results based on cluster-robust inference (var-cov estimator: CR1,
   approx t-test and confidence interval, df: residual method)
Code
round(i2_ml(mod_ma), 2)
    I2_Total  I2_FocalSpL     I2_RecNo I2_SubjectID    I2_Obs_ID 
       92.81         1.22        15.96         0.00        75.63 
Code
orchard_plot(mod_ma,
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code
# TODO
Code
# reduced model

mod_ma2 <- rma.mv(yi = SMD,
       V = VCV, 
       random = list(~1 | FocalSpL,
                     ~1 | RecNo,
                     ~1 | Obs_ID),
       test = "t",
       method = "REML", 
       sparse = TRUE,
       data = dat)

summary(mod_ma2)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-894.6453  1789.2906  1797.2906  1815.1302  1797.3537   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0111  0.1052     90     no  FocalSpL 
sigma^2.2  0.1450  0.3808    116     no     RecNo 
sigma^2.3  0.6872  0.8290    640     no    Obs_ID 

Model Results:

estimate      se    tval   df    pval   ci.lb   ci.ub      
  0.4139  0.0653  6.3401  639  <.0001  0.2857  0.5421  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(i2_ml(mod_ma2), 2)
   I2_Total I2_FocalSpL    I2_RecNo   I2_Obs_ID 
      92.81        1.22       15.96       75.63 
Code
# comparing two models
#anova(mod0, mod0r)

7 Meta-regression: uni-moderator

Code
mod_trt <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1 | FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mod = ~ Treat_mod - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod_trt)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-887.7086  1775.4172  1793.4172  1833.4856  1793.7056   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0074  0.0862     90     no  FocalSpL 
sigma^2.2  0.1563  0.3954    116     no     RecNo 
sigma^2.3  0.6907  0.8311    640     no    Obs_ID 

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 634) = 6.8763, p-val < .0001

Model Results:

              estimate      se    tval   df    pval    ci.lb   ci.ub      
Treat_modA      0.4019  0.0971  4.1381  634  <.0001   0.2112  0.5926  *** 
Treat_modAV     0.4723  0.1345  3.5128  634  0.0005   0.2083  0.7363  *** 
Treat_modAVO    0.5955  0.4487  1.3270  634  0.1850  -0.2857  1.4767      
Treat_modO      0.2577  0.2930  0.8797  634  0.3793  -0.3176  0.8330      
Treat_modOV     0.0185  0.3829  0.0482  634  0.9615  -0.7334  0.7703      
Treat_modV      0.4324  0.1043  4.1442  634  <.0001   0.2275  0.6373  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_trt)*100, 2)
   R2_marginal R2_conditional 
          0.64          19.68 
Code
orchard_plot(mod_trt, 
             mod = "Treat_mod",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code
# result table
all_models(mod_trt, mod = "Treat_mod")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
A 0.402 0.211 0.593 0.000 -1.423 2.227
AV 0.472 0.208 0.736 0.000 -1.362 2.307
AVO 0.595 -0.286 1.477 0.185 -1.422 2.613
O 0.258 -0.318 0.833 0.379 -1.646 2.162
OV 0.018 -0.733 0.770 0.962 -1.946 1.983
V 0.432 0.228 0.637 0.000 -1.394 2.259
A-AV 0.070 -0.248 0.389 0.664 -1.772 1.913
A-AVO 0.194 -0.706 1.093 0.673 -1.832 2.219
A-O -0.144 -0.748 0.460 0.639 -2.057 1.769
A-OV -0.383 -1.158 0.392 0.332 -2.357 1.590
A-V 0.031 -0.240 0.301 0.825 -1.805 1.866
AV-AVO 0.123 -0.797 1.043 0.793 -1.912 2.158
AV-O -0.215 -0.847 0.418 0.506 -2.137 1.708
AV-OV -0.454 -1.250 0.343 0.264 -2.436 1.528
AV-V -0.040 -0.340 0.260 0.794 -1.880 1.800
AVO-O -0.338 -1.390 0.714 0.529 -2.436 1.760
AVO-OV -0.577 -1.711 0.557 0.318 -2.717 1.563
AVO-V -0.163 -1.067 0.741 0.723 -2.191 1.865
O-OV -0.239 -1.186 0.707 0.620 -2.286 1.808
O-V 0.175 -0.434 0.784 0.574 -1.740 2.089
OV-V 0.414 -0.364 1.192 0.296 -1.561 2.389
Code
#homoscedastic model

VCVs <- vcalc(vi = dat_short$Vd,
             cluster = dat_short$SubjectID,
             rho = 0.5)


mod_trt_1 <- rma.mv(yi = SMD,
                   V = VCVs,
                   random = list(~1|FocalSpL,
                                 ~1 | RecNo,
                                 ~1 | Obs_ID),
                   mod = ~ Treat_mod - 1, 
                   test = "t",
                   method = "REML", 
                   sparse = TRUE,
                   data = dat_short)

summary(mod_trt_1)

Multivariate Meta-Analysis Model (k = 600; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-847.7054  1695.4108  1707.4108  1733.7623  1707.5532   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0087  0.0932     81     no  FocalSpL 
sigma^2.2  0.1700  0.4123    104     no     RecNo 
sigma^2.3  0.7138  0.8449    600     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 597) = 12.6006, p-val < .0001

Model Results:

             estimate      se    tval   df    pval   ci.lb   ci.ub      
Treat_modA     0.4055  0.0996  4.0715  597  <.0001  0.2099  0.6012  *** 
Treat_modV     0.4386  0.1069  4.1041  597  <.0001  0.2287  0.6485  *** 
Treat_modAV    0.4782  0.1379  3.4687  597  0.0006  0.2074  0.7490  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_trt_1)*100, 2)
   R2_marginal R2_conditional 
          0.08          20.09 
Code
orchard_plot(mod_trt_1, 
             mod = "Treat_mod",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code
# result table (get all the constrasts)
all_models(mod_trt_1, mod = "Treat_mod")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
A 0.406 0.210 0.601 0.000 -1.460 2.271
V 0.439 0.229 0.648 0.000 -1.429 2.306
AV 0.478 0.207 0.749 0.001 -1.397 2.353
A-V 0.033 -0.244 0.310 0.815 -1.843 1.909
A-AV 0.073 -0.254 0.399 0.662 -1.811 1.957
V-AV 0.040 -0.267 0.346 0.800 -1.841 1.920
Code
# modeling heteroscedasticity
mod_trt_2 <- rma.mv(yi = SMD, 
                V = VCVs, 
                random = list(~1|FocalSpL , 
                              ~1 | RecNo, 
                              ~ Treat_mod | Obs_ID), 
                mod = ~ Treat_mod - 1, 
                test = "t",
                struct = "DIAG",
                method = "REML", 
                sparse = TRUE,
                data = dat_short)

summary(mod_trt_2)

Multivariate Meta-Analysis Model (k = 600; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-838.2594  1676.5188  1692.5188  1727.6541  1692.7637   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0086  0.0927     81     no  FocalSpL 
sigma^2.2  0.1764  0.4200    104     no     RecNo 

outer factor: Obs_ID    (nlvls = 600)
inner factor: Treat_mod (nlvls = 3)

            estim    sqrt  k.lvl  fixed  level 
tau^2.1    0.7479  0.8648    302     no      A 
tau^2.2    0.8583  0.9264    190     no      V 
tau^2.3    0.3460  0.5882    108     no     AV 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 597) = 12.4198, p-val < .0001

Model Results:

             estimate      se    tval   df    pval   ci.lb   ci.ub      
Treat_modA     0.4100  0.1013  4.0495  597  <.0001  0.2112  0.6089  *** 
Treat_modV     0.4400  0.1108  3.9729  597  <.0001  0.2225  0.6576  *** 
Treat_modAV    0.4482  0.1229  3.6459  597  0.0003  0.2068  0.6897  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
orchard_plot(mod_trt_2, 
             mod = "Treat_mod",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code
# result table (get all the constrasts)
all_models(mod_trt_2, mod = "Treat_mod")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
A 0.410 0.211 0.609 0.000 -1.497 2.317
V 0.440 0.223 0.658 0.000 -1.578 2.458
AV 0.448 0.207 0.690 0.000 -1.003 1.900
A-V 0.033 -0.244 0.310 0.815 -1.843 1.909
A-AV 0.073 -0.254 0.399 0.662 -1.811 1.957
V-AV 0.040 -0.267 0.346 0.800 -1.841 1.920
Code
# comparison models
anova(mod_trt_1, mod_trt_2) # modeling heteroscedasticity is better 

        df       AIC       BIC      AICc    logLik     LRT   pval QE 
Full     8 1692.5188 1727.6541 1692.7637 -838.2594                NA 
Reduced  6 1707.4108 1733.7623 1707.5532 -847.7054 18.8921 <.0001 NA 
Code
# the effect of additions
# this is a part of sensitivity analysis

dat_short_add$Treat_add <- as.factor(paste(dat_short_add$Treat_mod , 
                                 dat_short_add$Add_on, 
                                 sep = "-"))

VCVs2 <- vcalc(vi = dat_short_add$Vd,
             cluster = dat_short_add$SubjectID,
             rho = 0.5)

mod_add <- rma.mv(yi = SMD, 
                V = VCVs2, 
                random = list(~1|FocalSpL , 
                              ~1 | RecNo, 
                              ~ 1 | Obs_ID), 
                mod = ~ Treat_add - 1,
                test = "t",
                struct = "DIAG",
                method = "REML", 
                sparse = TRUE,
                data = dat_short_add)

summary(mod_add) 

Multivariate Meta-Analysis Model (k = 298; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-422.9579   845.9158   859.9158   885.7008   860.3074   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0000     39     no  FocalSpL 
sigma^2.2  0.1432  0.3784     58     no     RecNo 
sigma^2.3  0.6666  0.8165    298     no    Obs_ID 

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 294) = 6.2038, p-val < .0001

Model Results:

                 estimate      se    tval   df    pval    ci.lb   ci.ub      
Treat_addAV-No     0.4270  0.1410  3.0289  294  0.0027   0.1495  0.7045   ** 
Treat_addAV-Yes    0.5716  0.2828  2.0210  294  0.0442   0.0150  1.1283    * 
Treat_addV-No      0.4262  0.1095  3.8912  294  0.0001   0.2106  0.6418  *** 
Treat_addV-Yes     0.2930  0.2091  1.4017  294  0.1621  -0.1184  0.7045      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
orchard_plot(mod_add, 
             mod = "Treat_add",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code
# result table (get all the constrasts)
all_models(mod_add, mod = "Treat_add")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
AV-No 0.427 0.150 0.704 0.003 -1.366 2.220
AV-Yes 0.572 0.015 1.128 0.044 -1.285 2.428
V-No 0.426 0.211 0.642 0.000 -1.358 2.210
V-Yes 0.293 -0.118 0.704 0.162 -1.525 2.111
AV-No-AV-Yes 0.145 -0.466 0.755 0.642 -1.729 2.018
AV-No-V-No -0.001 -0.324 0.322 0.996 -1.801 1.800
AV-No-V-Yes -0.134 -0.619 0.351 0.587 -1.970 1.702
AV-Yes-V-No -0.145 -0.734 0.443 0.627 -2.012 1.721
AV-Yes-V-Yes -0.279 -0.936 0.379 0.405 -2.168 1.611
V-No-V-Yes -0.133 -0.578 0.312 0.556 -1.959 1.693
Code
# testing the number of stimuli

mod_ord <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo, 
                             ~1 | Obs_ID), 
               mod = ~ Treat_No, 
               test = "t",
               method = "REML", 
               sparse = TRUE,
               data = dat)

summary(mod_ord)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-893.4118  1786.8236  1796.8236  1819.1153  1796.9185   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0117  0.1081     90     no  FocalSpL 
sigma^2.2  0.1469  0.3833    116     no     RecNo 
sigma^2.3  0.6875  0.8292    640     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 638) = 0.0952, p-val = 0.7577

Model Results:

          estimate      se    tval   df    pval    ci.lb   ci.ub    
intrcpt     0.3696  0.1603  2.3057  638  0.0214   0.0548  0.6843  * 
Treat_No    0.0366  0.1185  0.3086  638  0.7577  -0.1962  0.2694    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
bubble_plot(mod_ord,
             mod = "Treat_No",
             group = "RecNo",
             xlab = "The number of simuli",
             g = TRUE)

Code
#homoscedastic model
# Type of responses 
mod_type1 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1 | FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mod = ~ Type - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod_type1)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-886.7295  1773.4590  1785.4590  1812.1996  1785.5923   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0002     90     no  FocalSpL 
sigma^2.2  0.1322  0.3635    116     no     RecNo 
sigma^2.3  0.6902  0.8308    640     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 637) = 18.0966, p-val < .0001

Model Results:

                 estimate      se    tval   df    pval    ci.lb   ci.ub      
TypeBehaviour      0.4866  0.0677  7.1917  637  <.0001   0.3537  0.6194  *** 
TypeLifeHistory    0.3089  0.1176  2.6263  637  0.0088   0.0779  0.5398   ** 
TypePhysiology     0.0284  0.1285  0.2207  637  0.8254  -0.2240  0.2807      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_type1)*100, 2)
   R2_marginal R2_conditional 
          3.15          18.71 
Code
orchard_plot(mod_type1,
             mod = "Type",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code
# result table
all_models(mod_type1, mod = "Type")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Behaviour 0.487 0.354 0.619 0.000 -1.299 2.272
LifeHistory 0.309 0.078 0.540 0.009 -1.487 2.105
Physiology 0.028 -0.224 0.281 0.825 -1.770 1.827
Behaviour-LifeHistory -0.178 -0.418 0.062 0.146 -1.975 1.619
Behaviour-Physiology -0.458 -0.723 -0.193 0.001 -2.259 1.342
LifeHistory-Physiology -0.280 -0.586 0.025 0.072 -2.087 1.526
Code
# heteroscadasticity model
mod_type2 <- rma.mv(yi = SMD, 
               V = VCV, 
               mod = ~ Type - 1, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo, 
                             ~Type | Obs_ID), 
               struct = "DIAG",
               test = "t",
               method = "REML", 
               sparse = TRUE,
               data = dat)

summary(mod_type2)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-870.7865  1741.5731  1757.5731  1793.2272  1757.8024   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0274  0.1656     90     no  FocalSpL 
sigma^2.2  0.1296  0.3600    116     no     RecNo 

outer factor: Obs_ID (nlvls = 640)
inner factor: Type   (nlvls = 3)

            estim    sqrt  k.lvl  fixed        level 
tau^2.1    0.6577  0.8110    434     no    Behaviour 
tau^2.2    0.3176  0.5636    112     no  LifeHistory 
tau^2.3    1.2391  1.1131     94     no   Physiology 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 637) = 17.4552, p-val < .0001

Model Results:

                 estimate      se    tval   df    pval    ci.lb   ci.ub      
TypeBehaviour      0.5026  0.0720  6.9850  637  <.0001   0.3613  0.6439  *** 
TypeLifeHistory    0.3531  0.1047  3.3712  637  0.0008   0.1474  0.5587  *** 
TypePhysiology     0.0120  0.1549  0.0773  637  0.9384  -0.2922  0.3161      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
orchard_plot(mod_type2, 
             mod = "Type",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code
# result table
all_models(mod_type2, mod = "Type")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Behaviour 0.503 0.361 0.644 0.000 -1.275 2.281
LifeHistory 0.353 0.147 0.559 0.001 -1.015 1.721
Physiology 0.012 -0.292 0.316 0.938 -2.328 2.352
Behaviour-LifeHistory -0.178 -0.418 0.062 0.146 -1.975 1.619
Behaviour-Physiology -0.458 -0.723 -0.193 0.001 -2.259 1.342
LifeHistory-Physiology -0.280 -0.586 0.025 0.072 -2.087 1.526
Code
# heteroscadasticity model better than the homoscedasticity model
# note LifeHistory has lowest variation but this may be expected? 
# as it is less flexiable (e.g. the number of eggs?)
anova(mod_type1, mod_type2)

        df       AIC       BIC      AICc    logLik     LRT   pval QE 
Full     8 1757.5731 1793.2272 1757.8024 -870.7865                NA 
Reduced  6 1785.4590 1812.1996 1785.5923 -886.7295 31.8859 <.0001 NA 
Code
# treatment duration

dat$ln_duration <- log(dat$duration_days)

mod_dur <- rma.mv(yi = SMD,
               V = VCV,
               random = list(~1|FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ ln_duration,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod_dur)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-889.1410  1778.2821  1788.2821  1810.5738  1788.3770   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     90     no  FocalSpL 
sigma^2.2  0.1518  0.3897    116     no     RecNo 
sigma^2.3  0.6888  0.8300    640     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 638) = 8.4146, p-val = 0.0039

Model Results:

             estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt        0.3710  0.0634   5.8526  638  <.0001   0.2465   0.4955  *** 
ln_duration   -0.0455  0.0157  -2.9008  638  0.0039  -0.0763  -0.0147   ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_dur)*100, 2)
   R2_marginal R2_conditional 
          3.43          20.87 
Code
bubble_plot(mod_dur,
             mod = "ln_duration",
             group = "RecNo",
             xlab = "log(duration in days)",
             g = TRUE) +
    geom_point(data = dat,
    aes(x = ln_duration, y = SMD,
    color = Type,
    fill = Type,
    size = 1/sqrt(Vd)), alpha = 0.6) +
    scale_color_discrete() + 
    guides(color = "legend")

Code
#p + geom_point(aes(colour = Type))
#scale_colour_manual(values = c("red", "blue", "green"))

# extra  - looking at the interaction
mod_dur_int <- rma.mv(yi = SMD,
               V = VCV,
               random = list(~1|FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ ln_duration*Type,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

# no significant interaction
summary(mod_dur_int)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-880.1624  1760.3248  1778.3248  1818.3933  1778.6133   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     90     no  FocalSpL 
sigma^2.2  0.1427  0.3777    116     no     RecNo 
sigma^2.3  0.6882  0.8296    640     no    Obs_ID 

Test of Moderators (coefficients 2:6):
F(df1 = 5, df2 = 634) = 3.5544, p-val = 0.0035

Model Results:

                             estimate      se     tval   df    pval    ci.lb 
intrcpt                        0.4394  0.0717   6.1286  634  <.0001   0.2986 
ln_duration                   -0.0411  0.0173  -2.3802  634  0.0176  -0.0750 
TypeLifeHistory               -0.3341  0.3591  -0.9302  634  0.3526  -1.0393 
TypePhysiology                -0.4638  0.1637  -2.8326  634  0.0048  -0.7853 
ln_duration:TypeLifeHistory    0.0837  0.1012   0.8265  634  0.4088  -0.1151 
ln_duration:TypePhysiology     0.0438  0.0495   0.8849  634  0.3765  -0.0534 
                               ci.ub      
intrcpt                       0.5802  *** 
ln_duration                  -0.0072    * 
TypeLifeHistory               0.3712      
TypePhysiology               -0.1423   ** 
ln_duration:TypeLifeHistory   0.2825      
ln_duration:TypePhysiology    0.1409      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
res<-mod_results(mod_dur_int, mod = "ln_duration", group = "RecNo", by = "Type")

bubble_plot(res,
             mod = "ln_duration",
             group = "RecNo",
             xlab = "log(duration in days)",
             condition.nrow = 3,
             g = TRUE) 

Code
# sex

mod_sex1 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Sex - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod_sex1)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-890.5710  1781.1420  1793.1420  1819.8826  1793.2753   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0099  0.0993     90     no  FocalSpL 
sigma^2.2  0.1354  0.3680    116     no     RecNo 
sigma^2.3  0.6906  0.8310    640     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 637) = 15.0441, p-val < .0001

Model Results:

         estimate      se    tval   df    pval   ci.lb   ci.ub      
Sexboth    0.4706  0.0758  6.2110  637  <.0001  0.3218  0.6193  *** 
SexF       0.2579  0.1011  2.5501  637  0.0110  0.0593  0.4565    * 
SexM       0.4309  0.1325  3.2515  637  0.0012  0.1707  0.6912   ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_sex1)*100, 2)
   R2_marginal R2_conditional 
          0.95          18.17 
Code
orchard_plot(mod_sex1,
             mod = "Sex",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code
# result table
all_models(mod_sex1, mod = "Sex")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Both 0.471 0.322 0.619 0.000 -1.331 2.272
F 0.258 0.059 0.456 0.011 -1.548 2.064
M 0.431 0.171 0.691 0.001 -1.383 2.245
Both-F -0.213 -0.434 0.009 0.059 -2.022 1.596
Both-M -0.040 -0.323 0.244 0.784 -1.857 1.778
F-M 0.173 -0.098 0.444 0.211 -1.643 1.989
Code
# heteroscedasticity

mod_sex2 <- rma.mv(yi = SMD, 
                V = VCV, 
                mod = ~ Sex - 1, 
                random = list(~1|FocalSpL , 
                              ~1 | RecNo, 
                              ~Sex | Obs_ID), 
                struct = "DIAG",
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dat)

summary(mod_sex2)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-867.9399  1735.8799  1751.8799  1787.5341  1752.1092   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0007     90     no  FocalSpL 
sigma^2.2  0.1454  0.3813    116     no     RecNo 

outer factor: Obs_ID (nlvls = 640)
inner factor: Sex    (nlvls = 3)

            estim    sqrt  k.lvl  fixed  level 
tau^2.1    0.8865  0.9415    398     no   both 
tau^2.2    0.2883  0.5369    158     no      F 
tau^2.3    0.6740  0.8210     84     no      M 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 637) = 14.8707, p-val < .0001

Model Results:

         estimate      se    tval   df    pval   ci.lb   ci.ub      
Sexboth    0.4729  0.0781  6.0557  637  <.0001  0.3196  0.6263  *** 
SexF       0.2429  0.0834  2.9127  637  0.0037  0.0791  0.4066   ** 
SexM       0.4268  0.1278  3.3403  637  0.0009  0.1759  0.6777  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
orchard_plot(mod_sex2,
             mod = "Sex",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code
# result table
all_models(mod_sex2, mod = "Sex", type = "hetero")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Both 0.473 0.320 0.626 0.000 -1.528 2.474
F 0.243 0.079 0.407 0.004 -1.061 1.546
M 0.427 0.176 0.678 0.001 -1.368 2.222
Both-F -0.213 -0.434 0.009 0.059 -2.022 1.596
Both-M -0.040 -0.323 0.244 0.784 -1.857 1.778
F-M 0.173 -0.098 0.444 0.211 -1.643 1.989
Code
anova(mod_sex1, mod_sex2)

        df       AIC       BIC      AICc    logLik     LRT   pval QE 
Full     8 1751.8799 1787.5341 1752.1092 -867.9399                NA 
Reduced  6 1793.1420 1819.8826 1793.2753 -890.5710 45.2621 <.0001 NA 
Code
# Category of responses


mod_cat <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo, 
                             ~1 | Obs_ID), 
               mod = ~ Category - 1, 
               test = "t",
               method = "REML", 
               sparse = TRUE,
               data = dat)

summary(mod_cat)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-879.5908  1759.1816  1781.1816  1830.1194  1781.6074   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     90     no  FocalSpL 
sigma^2.2  0.1389  0.3727    116     no     RecNo 
sigma^2.3  0.6912  0.8314    640     no    Obs_ID 

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 632) = 7.2433, p-val < .0001

Model Results:

                          estimate      se    tval   df    pval    ci.lb 
CategoryAntiPredator        0.4772  0.0924  5.1669  632  <.0001   0.2959 
CategoryCondition           0.0027  0.1394  0.0191  632  0.9848  -0.2710 
CategoryCostlyBehaviours    0.4567  0.1784  2.5596  632  0.0107   0.1063 
CategoryHormones            0.1717  0.2914  0.5891  632  0.5560  -0.4006 
CategoryIntake              0.7554  0.1623  4.6539  632  <.0001   0.4366 
CategoryParentalCare        0.4103  0.1076  3.8133  632  0.0002   0.1990 
CategoryPhenology           0.1927  0.1799  1.0711  632  0.2845  -0.1606 
Categoryreproduction        0.3366  0.1287  2.6159  632  0.0091   0.0839 
                           ci.ub      
CategoryAntiPredator      0.6586  *** 
CategoryCondition         0.2763      
CategoryCostlyBehaviours  0.8071    * 
CategoryHormones          0.7439      
CategoryIntake            1.0741  *** 
CategoryParentalCare      0.6216  *** 
CategoryPhenology         0.5459      
Categoryreproduction      0.5892   ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_cat)*100, 2)
   R2_marginal R2_conditional 
          4.12          20.17 
Code
orchard_plot(mod_cat, 
             mod = "Category",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             angle = 45,
             branch.size = 3)

Code
# result table
all_models(mod_cat, mod = "Category")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
AntiPredator 0.477 0.296 0.659 0.000 -1.321 2.276
Condition 0.003 -0.271 0.276 0.985 -1.807 1.813
CostlyBehaviours 0.457 0.106 0.807 0.011 -1.366 2.280
Hormones 0.172 -0.401 0.744 0.556 -1.707 2.050
Intake 0.755 0.437 1.074 0.000 -1.062 2.573
ParentalCare 0.410 0.199 0.622 0.000 -1.391 2.212
Phenology 0.193 -0.161 0.546 0.285 -1.631 2.016
Reproduction 0.337 0.084 0.589 0.009 -1.470 2.143
AntiPredator-Condition -0.475 -0.788 -0.161 0.003 -2.291 1.342
AntiPredator-CostlyBehaviours -0.021 -0.389 0.348 0.913 -1.847 1.806
AntiPredator-Hormones -0.306 -0.901 0.290 0.314 -2.191 1.580
AntiPredator-Intake 0.278 -0.063 0.619 0.110 -1.543 2.100
AntiPredator-ParentalCare -0.067 -0.329 0.195 0.616 -1.875 1.741
AntiPredator-Phenology -0.285 -0.662 0.093 0.140 -2.113 1.544
AntiPredator-Reproduction -0.141 -0.428 0.146 0.336 -1.953 1.671
Condition-CostlyBehaviours 0.454 0.018 0.891 0.042 -1.388 2.296
Condition-Hormones 0.169 -0.445 0.783 0.589 -1.723 2.061
Condition-Intake 0.753 0.353 1.152 0.000 -1.081 2.586
Condition-ParentalCare 0.408 0.086 0.730 0.013 -1.410 2.226
Condition-Phenology 0.190 -0.228 0.608 0.373 -1.647 2.027
Condition-Reproduction 0.334 -0.004 0.671 0.053 -1.487 2.155
CostlyBehaviours-Hormones -0.285 -0.952 0.382 0.402 -2.195 1.625
CostlyBehaviours-Intake 0.299 -0.160 0.757 0.202 -1.548 2.146
CostlyBehaviours-ParentalCare -0.046 -0.448 0.356 0.821 -1.880 1.787
CostlyBehaviours-Phenology -0.264 -0.753 0.225 0.289 -2.119 1.591
CostlyBehaviours-Reproduction -0.120 -0.543 0.303 0.578 -1.959 1.718
Hormones-Intake 0.584 -0.068 1.235 0.079 -1.320 2.488
Hormones-ParentalCare 0.239 -0.362 0.839 0.435 -1.649 2.126
Hormones-Phenology 0.021 -0.625 0.667 0.949 -1.881 1.923
Hormones-Reproduction 0.165 -0.447 0.777 0.597 -1.726 2.056
Intake-ParentalCare -0.345 -0.720 0.030 0.071 -2.173 1.483
Intake-Phenology -0.563 -1.031 -0.094 0.019 -2.412 1.287
Intake-Reproduction -0.419 -0.817 -0.020 0.039 -2.252 1.414
ParentalCare-Phenology -0.218 -0.611 0.176 0.278 -2.050 1.614
ParentalCare-Reproduction -0.074 -0.374 0.227 0.630 -1.888 1.740
Phenology-Reproduction 0.144 -0.220 0.508 0.438 -1.682 1.970
Code
# Redactor guild
# Fish only N = 1 so remove and also remove NS (not specified)

dat %>% filter(PredGuild != "Fish" & PredGuild != "NS") ->
  data_preg

VCVs3 <- vcalc(vi = data_preg$Vd,
             cluster = data_preg$SubjectID,
             rho = 0.5)

mod_preg <- rma.mv(yi = SMD, 
               V = VCVs3, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ PredGuild - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = data_preg)

summary(mod_preg)

Multivariate Meta-Analysis Model (k = 620; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-857.4707  1714.9413  1728.9413  1759.9040  1729.1255   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0007  0.0271     85     no  FocalSpL 
sigma^2.2  0.1517  0.3895    110     no     RecNo 
sigma^2.3  0.6839  0.8270    620     no    Obs_ID 

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 616) = 10.0849, p-val < .0001

Model Results:

                  estimate      se    tval   df    pval    ci.lb   ci.ub      
PredGuildBird       0.4520  0.0732  6.1774  616  <.0001   0.3083  0.5957  *** 
PredGuildMammal     0.2681  0.1279  2.0957  616  0.0365   0.0169  0.5193    * 
PredGuildmixed      0.1120  0.2485  0.4508  616  0.6523  -0.3760  0.6001      
PredGuildReptile    0.3220  0.2800  1.1501  616  0.2506  -0.2279  0.8719      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_preg)*100, 2)
   R2_marginal R2_conditional 
          1.70          19.62 
Code
# R2_marginal R2_conditional 
#           1.70          19.62 

orchard_plot(mod_preg, 
             mod = "PredGuild",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code
# result table
all_models(mod_preg, mod = "PredGuild")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Bird 0.452 0.308 0.596 0.000 -1.350 2.254
Mammal 0.268 0.017 0.519 0.037 -1.545 2.082
Mixed 0.112 -0.376 0.600 0.652 -1.749 1.973
Reptile 0.322 -0.228 0.872 0.251 -1.556 2.200
Bird-Mammal -0.184 -0.451 0.083 0.176 -2.000 1.632
Bird-Mixed -0.340 -0.847 0.167 0.189 -2.206 1.526
Bird-Reptile -0.130 -0.695 0.435 0.651 -2.013 1.753
Mammal-Mixed -0.156 -0.705 0.392 0.577 -2.034 1.722
Mammal-Reptile 0.054 -0.550 0.657 0.861 -1.841 1.949
Mixed-Reptile 0.210 -0.525 0.945 0.575 -1.731 2.151
Code
# Setting

mod_set <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Setting - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod_set)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-891.0585  1782.1169  1794.1169  1820.8575  1794.2502   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0172  0.1311     90     no  FocalSpL 
sigma^2.2  0.1378  0.3712    116     no     RecNo 
sigma^2.3  0.6868  0.8288    640     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 637) = 14.2658, p-val < .0001

Model Results:

                     estimate      se     tval   df    pval    ci.lb   ci.ub 
SettingField           0.4245  0.0721   5.8846  637  <.0001   0.2829  0.5662 
SettingLab             0.4764  0.1575   3.0247  637  0.0026   0.1671  0.7857 
SettingSemi-natural   -0.2054  0.4204  -0.4886  637  0.6253  -1.0308  0.6201 
                         
SettingField         *** 
SettingLab            ** 
SettingSemi-natural      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_set)*100, 2)
   R2_marginal R2_conditional 
          0.70          18.98 
Code
# R2_marginal R2_conditional 
#           0.70          18.98 

orchard_plot(mod_set, 
             mod = "Setting",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code
# result table
all_models(mod_set, mod = "Setting")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Field 0.425 0.283 0.566 0.000 -1.383 2.232
Lab 0.476 0.167 0.786 0.003 -1.352 2.304
Semi-natural -0.205 -1.031 0.620 0.625 -2.187 1.776
Field-Lab 0.052 -0.283 0.387 0.761 -1.781 1.884
Field-Semi-natural -0.630 -1.463 0.204 0.138 -2.615 1.355
Lab-Semi-natural -0.682 -1.559 0.195 0.127 -2.686 1.322
Code
# Season
mod_4 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Season - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)


summary(mod_4)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-892.1745  1784.3489  1794.3489  1816.6406  1794.4439   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0187  0.1368     90     no  FocalSpL 
sigma^2.2  0.1390  0.3728    116     no     RecNo 
sigma^2.3  0.6875  0.8291    640     no    Obs_ID 

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 638) = 21.2697, p-val < .0001

Model Results:

                    estimate      se    tval   df    pval   ci.lb   ci.ub      
Seasonbreeding        0.3589  0.0783  4.5831  638  <.0001  0.2051  0.5127  *** 
Seasonnon-breeding    0.5596  0.1123  4.9845  638  <.0001  0.3392  0.7801  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_4)*100, 2)
   R2_marginal R2_conditional 
          0.95          19.44 
Code
orchard_plot(mod_4,
             mod = "Season",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code
# result table
all_models(mod_4, mod = "Season")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Breeding 0.359 0.205 0.513 0.000 -1.453 2.171
Non-breeding 0.560 0.339 0.780 0.000 -1.259 2.378
Breeding-Non-breeding 0.201 -0.058 0.460 0.128 -1.623 2.025
Code
# Design
mod_d <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Design - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod_d)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-893.3891  1786.7783  1796.7783  1819.0700  1796.8732   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0103  0.1016     90     no  FocalSpL 
sigma^2.2  0.1481  0.3848    116     no     RecNo 
sigma^2.3  0.6879  0.8294    640     no    Obs_ID 

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 638) = 19.9969, p-val < .0001

Model Results:

              estimate      se    tval   df    pval   ci.lb   ci.ub      
Designamong     0.4061  0.0874  4.6483  638  <.0001  0.2345  0.5777  *** 
Designwithin    0.4197  0.0839  5.0006  638  <.0001  0.2549  0.5844  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_d)*100, 2)
   R2_marginal R2_conditional 
          0.01          18.72 
Code
 # R2_marginal R2_conditional 
 #          0.01          18.72 

orchard_plot(mod_d,
             mod = "Design",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code
# result table
all_models(mod_d, mod = "Design")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Among 0.406 0.235 0.578 0.000 -1.409 2.221
Within 0.420 0.255 0.584 0.000 -1.394 2.234
Among-Within 0.014 -0.204 0.231 0.903 -1.806 1.833
Code
# Response period
mod_resp <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ ResponsePeriod - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod_resp)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-892.5749  1785.1498  1795.1498  1817.4415  1795.2448   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0138  0.1176     90     no  FocalSpL 
sigma^2.2  0.1435  0.3788    116     no     RecNo 
sigma^2.3  0.6883  0.8296    640     no    Obs_ID 

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 638) = 20.8414, p-val < .0001

Model Results:

                      estimate      se    tval   df    pval   ci.lb   ci.ub 
ResponsePeriodafter     0.3500  0.0851  4.1108  638  <.0001  0.1828  0.5171 
ResponsePeriodduring    0.4938  0.0903  5.4708  638  <.0001  0.3166  0.6711 
                          
ResponsePeriodafter   *** 
ResponsePeriodduring  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_resp)*100, 2)
   R2_marginal R2_conditional 
          0.61          19.10 
Code
#   R2_marginal R2_conditional 
#          0.61          19.10 


orchard_plot(mod_resp,
             mod = "ResponsePeriod",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code
# result table
all_models(mod_resp, mod = "ResponsePeriod")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
After 0.350 0.183 0.517 0.000 -1.463 2.163
During 0.494 0.317 0.671 0.000 -1.321 2.308
After-During 0.144 -0.083 0.371 0.214 -1.676 1.964
Code
# control type
# Mix only N = 1 so remove 

dat %>% filter(ControlType != "mix") ->
  data_cont

VCVs4 <- vcalc(vi = data_cont$Vd,
             cluster = data_cont$SubjectID,
             rho = 0.5)

mod_cont <- rma.mv(yi = SMD, 
               V = VCVs4, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ ControlType - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = data_cont)
summary(mod_cont)

Multivariate Meta-Analysis Model (k = 639; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-891.3557  1782.7115  1794.7115  1821.4427  1794.8450   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0080  0.0893     90     no  FocalSpL 
sigma^2.2  0.1504  0.3879    115     no     RecNo 
sigma^2.3  0.6901  0.8307    639     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 636) = 13.3306, p-val < .0001

Model Results:

                        estimate      se    tval   df    pval   ci.lb   ci.ub 
ControlTypeBlank          0.4424  0.0947  4.6709  636  <.0001  0.2564  0.6283 
ControlTypedisturbance    0.4170  0.1452  2.8710  636  0.0042  0.1318  0.7022 
ControlTypeNonPred        0.3895  0.0810  4.8106  636  <.0001  0.2305  0.5485 
                            
ControlTypeBlank        *** 
ControlTypedisturbance   ** 
ControlTypeNonPred      *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_cont)*100, 2)
   R2_marginal R2_conditional 
          0.07          18.73 
Code
# R2_marginal R2_conditional 
#           0.07          18.73 

orchard_plot(mod_cont,
             mod = "ControlType",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code
# result table
all_models(mod_cont, mod = "ControlType")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Blank 0.442 0.256 0.628 0.000 -1.376 2.261
Disturbance 0.417 0.132 0.702 0.004 -1.414 2.248
NonPred 0.389 0.231 0.548 0.000 -1.426 2.205
Blank-Disturbance -0.025 -0.346 0.295 0.876 -1.862 1.812
Blank-NonPred -0.053 -0.260 0.155 0.617 -1.874 1.768
Disturbance-NonPred -0.027 -0.346 0.291 0.865 -1.864 1.809
Code
# age
#A = Adult, N = Nestling, J = Juveniles (including fledglings and first year birds), E = Egg

# "" only N = 2 (study = 1) so remove 

dat %>% filter(Age != "") ->
  dat_age

VCVs5 <- vcalc(vi = dat_age$Vd,
             cluster = dat_age$SubjectID,
             rho = 0.5)


# age
mod_age <- rma.mv(yi = SMD, 
               V = VCVs5, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Age - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat_age)
summary(mod_age)

Multivariate Meta-Analysis Model (k = 638; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-888.1737  1776.3474  1790.3474  1821.5117  1790.5263   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0084  0.0917     90     no  FocalSpL 
sigma^2.2  0.1498  0.3870    116     no     RecNo 
sigma^2.3  0.6914  0.8315    638     no    Obs_ID 

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 634) = 9.9438, p-val < .0001

Model Results:

      estimate      se    tval   df    pval    ci.lb   ci.ub      
AgeA    0.4206  0.0677  6.2136  634  <.0001   0.2877  0.5536  *** 
AgeE    0.3228  0.3452  0.9353  634  0.3500  -0.3550  1.0006      
AgeJ    0.2881  0.3484  0.8269  634  0.4086  -0.3961  0.9724      
AgeN    0.3265  0.1911  1.7087  634  0.0880  -0.0487  0.7018    . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_age)*100, 2)
   R2_marginal R2_conditional 
          0.15          18.74 
Code
   # R2_marginal R2_conditional 
   #        0.15          18.74 

orchard_plot(mod_age,
             mod = "Age",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code
# result table
all_models(mod_age, mod = "Age")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
A 0.421 0.288 0.554 0.000 -1.394 2.235
E 0.323 -0.355 1.001 0.350 -1.610 2.256
J 0.288 -0.396 0.972 0.409 -1.647 2.223
N 0.327 -0.049 0.702 0.088 -1.522 2.175
A-E -0.098 -0.782 0.587 0.779 -2.033 1.837
A-J -0.133 -0.823 0.558 0.706 -2.070 1.805
A-N -0.094 -0.472 0.283 0.625 -1.943 1.755
E-J -0.035 -0.997 0.927 0.944 -2.085 2.015
E-N 0.004 -0.766 0.774 0.992 -1.963 1.971
J-N 0.038 -0.738 0.815 0.923 -1.931 2.008
Code
# type of predator

dat$PredTo[dat$PredTo == ""] <- NA



mod_predt <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ PredTo - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod_predt)

Multivariate Meta-Analysis Model (k = 621; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-861.2315  1722.4629  1734.4629  1761.0219  1734.6004   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0003     85     no  FocalSpL 
sigma^2.2  0.1572  0.3965    110     no     RecNo 
sigma^2.3  0.6805  0.8249    621     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 618) = 12.3559, p-val < .0001

Model Results:

         estimate      se    tval   df    pval    ci.lb   ci.ub      
PredToA    0.4284  0.0814  5.2628  618  <.0001   0.2686  0.5883  *** 
PredToB    0.3772  0.3065  1.2306  618  0.2190  -0.2247  0.9791      
PredToN    0.3337  0.0915  3.6466  618  0.0003   0.1540  0.5135  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_predt)*100, 2)
   R2_marginal R2_conditional 
          0.25          18.97 
Code
# R2_marginal R2_conditional 
#           0.25          18.97 


orchard_plot(mod_predt,
             mod = "PredTo",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code
# result table
all_models(mod_predt, mod = "PredTo")
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
A 0.428 0.269 0.588 0.000 -1.376 2.233
B 0.377 -0.225 0.979 0.219 -1.518 2.273
N 0.334 0.154 0.513 0.000 -1.473 2.140
A-B -0.051 -0.674 0.572 0.872 -1.953 1.851
A-N -0.095 -0.314 0.125 0.397 -1.905 1.716
B-N -0.043 -0.672 0.585 0.892 -1.947 1.861

8 Meta-regression: multi-moderator

Code
#######################
# Mulit-variable models
#######################

dat_short$sln_duration <- scale(log(dat_short$duration_days))

mod_full <- rma.mv(yi = SMD, 
               V = VCVs, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ sln_duration*Type +
                         sln_duration*Treat_mod +
                         Sex,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat_short)
summary(mod_full)

Multivariate Meta-Analysis Model (k = 600; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-826.1989  1652.3977  1682.3977  1748.0486  1683.2369   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     81     no  FocalSpL 
sigma^2.2  0.1540  0.3925    104     no     RecNo 
sigma^2.3  0.7169  0.8467    600     no    Obs_ID 

Test of Moderators (coefficients 2:12):
F(df1 = 11, df2 = 588) = 2.0578, p-val = 0.0216

Model Results:

                              estimate      se     tval   df    pval    ci.lb 
intrcpt                         0.4853  0.1076   4.5122  588  <.0001   0.2741 
sln_duration                   -0.1333  0.0896  -1.4886  588  0.1371  -0.3093 
TypeLifeHistory                -0.0267  0.4683  -0.0569  588  0.9546  -0.9463 
TypePhysiology                 -0.5076  0.1840  -2.7580  588  0.0060  -0.8691 
Treat_modV                      0.0477  0.1461   0.3264  588  0.7442  -0.2393 
Treat_modAV                     0.1566  0.1656   0.9455  588  0.3448  -0.1687 
SexF                           -0.2748  0.1295  -2.1219  588  0.0343  -0.5291 
SexM                           -0.0870  0.1531  -0.5679  588  0.5703  -0.3877 
sln_duration:TypeLifeHistory    0.1502  0.4474   0.3358  588  0.7372  -0.7284 
sln_duration:TypePhysiology     0.2185  0.1995   1.0952  588  0.2739  -0.1733 
sln_duration:Treat_modV        -0.0308  0.1418  -0.2168  588  0.8284  -0.3093 
sln_duration:Treat_modAV       -0.1650  0.1804  -0.9150  588  0.3606  -0.5193 
                                ci.ub      
intrcpt                        0.6966  *** 
sln_duration                   0.0426      
TypeLifeHistory                0.8930      
TypePhysiology                -0.1461   ** 
Treat_modV                     0.3347      
Treat_modAV                    0.4818      
SexF                          -0.0205    * 
SexM                           0.2138      
sln_duration:TypeLifeHistory   1.0289      
sln_duration:TypePhysiology    0.6103      
sln_duration:Treat_modV        0.2478      
sln_duration:Treat_modAV       0.1892      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_full)*100, 2)
   R2_marginal R2_conditional 
          7.06          23.50 
Code
orchard_plot(mod_full,
             mod = "Type",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code
orchard_plot(mod_full,
             mod = "Treat_mod",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code
# interaction plot of duration with trait types
int_type <- mod_results(mod_full, mod = "sln_duration", group = "RecNo", weights = "prop",
                                   by = "Type")

bubble_plot(int_type, group = "RecNo", mod = "sln_duration", xlab = "ln(duration in days) (z-transformed)",
                     legend.pos = "top.left", condition.nrow = 3, g = TRUE)

Code
# interaction plot of duration with trait types
int_trt <- mod_results(mod_full, mod = "sln_duration", group = "RecNo", weights = "prop",
                        by = "Treat_mod")

bubble_plot(int_trt, group = "RecNo", mod = "sln_duration", xlab = "ln(duration in days)",
            legend.pos = "top.left", condition.nrow = 3, g = TRUE)

Code
# mulit-model selection
candidates <- dredge(mod_full, trace = 2)

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=                                                                     |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |=======                                                               |   9%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |==================================================                    |  72%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |===================================================================== |  98%
Code
# displays delta AICc <2
candidates_aic2 <- subset(candidates, delta < 5) 
# model averaging
mr_averaged_aic2 <- summary(model.avg(candidates, delta < 5)) 

# relative importance of each predictor for all the models
importance <- sw(candidates)

9 Publication bias

Code
funnel(mod_ma2, 
       yaxis="seinv",
       # = "rstudent",
       xlab = "Standarized residuals",
       ylab = "Precision (inverse of SE)",
       ylim = c(0.001, 16),
       xlim = c(-10,15),
       digits=c(0,1)
       )

Code
funnel(mod_full, 
       yaxis="seinv",
       #type = "rstudent",
       xlab = "Standarized residuals",
       ylab = "Precision (inverse of SE)",
       ylim = c(0.001, 1.5),
       xlim = c(-5,10),
       digits=c(0,1)
       )

Code
# Egger

dat$effectN <- (dat$Ncontrol * dat$NTreat) / (dat$Ncontrol + dat$NTreat)
dat$sqeffectN <- sqrt(dat$effectN)

mod_mae <- rma.mv(yi = SMD,
                V = VCV,
                mods = ~ sqeffectN,
                random = list(
                  ~1 | FocalSpL,
                  ~1 | RecNo,
                  ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dat)

summary(mod_mae)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-893.0018  1786.0037  1796.0037  1818.2954  1796.0986   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0003     90     no  FocalSpL 
sigma^2.2  0.1483  0.3850    116     no     RecNo 
sigma^2.3  0.6929  0.8324    640     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 638) = 1.0650, p-val = 0.3025

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub      
intrcpt      0.5044  0.1194   4.2234  638  <.0001   0.2699  0.7390  *** 
sqeffectN   -0.0291  0.0282  -1.0320  638  0.3025  -0.0844  0.0262      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
bubble_plot(mod_mae,
            mod = "sqeffectN",
            group = "RecNo",
            xlab = "Effective N",
            g = TRUE)

Code
# decline effect
mod_mad <- rma.mv(yi = SMD,
                V = VCV,
                mods = ~ Year,
                random = list(
                  ~1 | FocalSpL,
                  ~1 | RecNo,
                  ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dat)

summary(mod_mad)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-892.5966  1785.1932  1795.1932  1817.4849  1795.2881   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0103  0.1016     90     no  FocalSpL 
sigma^2.2  0.1450  0.3808    116     no     RecNo 
sigma^2.3  0.6886  0.8298    640     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 638) = 1.5062, p-val = 0.2202

Model Results:

         estimate       se     tval   df    pval     ci.lb    ci.ub    
intrcpt   23.4292  18.7537   1.2493  638  0.2120  -13.3973  60.2556    
Year      -0.0114   0.0093  -1.2273  638  0.2202   -0.0297   0.0069    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
bubble_plot(mod_mad,
            mod = "Year",
            group = "RecNo",
            xlab = "Publication year",
            g = TRUE)

Code
# full model
dat_short$effectN <- (dat_short$Ncontrol * dat_short$NTreat) / (dat_short$Ncontrol + dat_short$NTreat)
dat_short$sqeffectN <- sqrt(dat_short$effectN)

mod_fulle <- rma.mv(yi = SMD, 
                   V = VCVs, 
                   random = list(~1|FocalSpL , 
                                 ~1 | RecNo,
                                 ~1 | Obs_ID),
                   mods =  ~ sqeffectN +
                     Year +
                     sln_duration*Type +
                     sln_duration*Treat_mod +
                     Sex,
                   test = "t",
                   method = "REML",
                   sparse = TRUE,
                   data = dat_short)
summary(mod_fulle)

Multivariate Meta-Analysis Model (k = 600; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-823.4875  1646.9750  1680.9750  1755.3215  1682.0525   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0002     81     no  FocalSpL 
sigma^2.2  0.1587  0.3983    104     no     RecNo 
sigma^2.3  0.7183  0.8475    600     no    Obs_ID 

Test of Moderators (coefficients 2:14):
F(df1 = 13, df2 = 586) = 1.7639, p-val = 0.0454

Model Results:

                              estimate       se     tval   df    pval     ci.lb 
intrcpt                        14.6850  22.7802   0.6446  586  0.5194  -30.0557 
sqeffectN                      -0.0007   0.0328  -0.0216  586  0.9827   -0.0651 
Year                           -0.0070   0.0113  -0.6221  586  0.5341   -0.0293 
sln_duration                   -0.1252   0.0936  -1.3370  586  0.1817   -0.3090 
TypeLifeHistory                -0.0173   0.4714  -0.0366  586  0.9708   -0.9431 
TypePhysiology                 -0.5322   0.1894  -2.8094  586  0.0051   -0.9043 
Treat_modV                      0.0368   0.1478   0.2489  586  0.8035   -0.2535 
Treat_modAV                     0.1375   0.1690   0.8136  586  0.4162   -0.1944 
SexF                           -0.2740   0.1300  -2.1076  586  0.0355   -0.5293 
SexM                           -0.0953   0.1543  -0.6178  586  0.5369   -0.3983 
sln_duration:TypeLifeHistory    0.1450   0.4510   0.3216  586  0.7479   -0.7407 
sln_duration:TypePhysiology     0.2438   0.2055   1.1863  586  0.2360   -0.1599 
sln_duration:Treat_modV        -0.0261   0.1440  -0.1813  586  0.8562   -0.3088 
sln_duration:Treat_modAV       -0.1497   0.1834  -0.8161  586  0.4148   -0.5099 
                                ci.ub     
intrcpt                       59.4257     
sqeffectN                      0.0637     
Year                           0.0152     
sln_duration                   0.0587     
TypeLifeHistory                0.9086     
TypePhysiology                -0.1602  ** 
Treat_modV                     0.3271     
Treat_modAV                    0.4693     
SexF                          -0.0187   * 
SexM                           0.2077     
sln_duration:TypeLifeHistory   1.0307     
sln_duration:TypePhysiology    0.6475     
sln_duration:Treat_modV        0.2566     
sln_duration:Treat_modAV       0.2105     

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# small-study
bubble_plot(mod_fulle, group = "RecNo", mod = "sqeffectN", xlab = "ln(duration in days)",
            legend.pos = "top.left", condition.nrow = 3, g = TRUE)

Code
# decline
bubble_plot(mod_fulle, group = "RecNo", mod = "Year", xlab = "ln(duration in days)",
            legend.pos = "top.left", condition.nrow = 3, g = TRUE)

11 R Session Informtion

Code
# pander for making it look nicer
sessionInfo() %>% pander()

R version 4.3.1 (2023-06-16)

Platform: aarch64-apple-darwin20 (64-bit)

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: pacman(v.0.5.1), pander(v.0.6.5), kableExtra(v.1.3.4), MuMIn(v.1.47.5), emmeans(v.1.8.7), clubSandwich(v.0.5.10), ape(v.5.7-1), easyalluvial(v.0.3.1), ggalluvial(v.0.12.5), alluvial(v.0.1-2), patchwork(v.1.1.3), metafor(v.4.2-0), numDeriv(v.2016.8-1.1), metadat(v.1.2-0), orchaRd(v.2.0), lme4(v.1.1-34), Matrix(v.1.6-1), here(v.1.0.1), lubridate(v.1.9.2), forcats(v.1.0.0), stringr(v.1.5.0), dplyr(v.1.1.2), purrr(v.1.0.2), readr(v.2.1.4), tidyr(v.1.3.0), tibble(v.3.2.1), ggplot2(v.3.4.3) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): gridExtra(v.2.3), sandwich(v.3.0-2), rlang(v.1.1.1), magrittr(v.2.0.3), ggridges(v.0.5.4), compiler(v.4.3.1), mgcv(v.1.8-42), systemfonts(v.1.0.4), vctrs(v.0.6.3), rvest(v.1.0.3), pkgconfig(v.2.0.3), fastmap(v.1.1.1), labeling(v.0.4.2), utf8(v.1.2.3), rmarkdown(v.2.24), prodlim(v.2023.03.31), tzdb(v.0.4.0), ggbeeswarm(v.0.7.2), nloptr(v.2.0.3), xfun(v.0.40), jsonlite(v.1.8.7), recipes(v.1.0.7), highr(v.0.10), parallel(v.4.3.1), R6(v.2.5.1), stringi(v.1.7.12), RColorBrewer(v.1.1-3), parallelly(v.1.36.0), boot(v.1.3-28.1), rpart(v.4.1.19), estimability(v.1.4.1), Rcpp(v.1.0.11), knitr(v.1.43), future.apply(v.1.11.0), zoo(v.1.8-12), splines(v.4.3.1), nnet(v.7.3-19), timechange(v.0.2.0), tidyselect(v.1.2.0), rstudioapi(v.0.15.0), yaml(v.2.3.7), timeDate(v.4022.108), codetools(v.0.2-19), listenv(v.0.9.0), lattice(v.0.21-8), withr(v.2.5.0), coda(v.0.19-4), evaluate(v.0.21), future(v.1.33.0), survival(v.3.5-5), xml2(v.1.3.5), pillar(v.1.9.0), stats4(v.4.3.1), generics(v.0.1.3), rprojroot(v.2.0.3), mathjaxr(v.1.6-0), hms(v.1.1.3), munsell(v.0.5.0), scales(v.1.2.1), minqa(v.1.2.5), globals(v.0.16.2), xtable(v.1.8-4), class(v.7.3-22), glue(v.1.6.2), tools(v.4.3.1), data.table(v.1.14.8), webshot(v.0.5.5), gower(v.1.0.1), mvtnorm(v.1.2-2), grid(v.4.3.1), ipred(v.0.9-14), colorspace(v.2.1-0), nlme(v.3.1-162), beeswarm(v.0.4.0), vipor(v.0.4.5), latex2exp(v.0.9.6), cli(v.3.6.1), fansi(v.1.0.4), viridisLite(v.0.4.2), svglite(v.2.1.1), lava(v.1.7.2.1), gtable(v.0.3.3), digest(v.0.6.33), progressr(v.0.14.0), farver(v.2.1.1), htmlwidgets(v.1.6.2), htmltools(v.0.5.6), lifecycle(v.1.0.3), hardhat(v.1.3.0), httr(v.1.4.7) and MASS(v.7.3-60)